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ABSTRACT 
In  an  input-output  econometric  model  we  are  often  concerned  with 
solving  the  system  of  n  equations  (I  -  A)  x  =  y,  repeatedly  for  various 
changes  in  the  elements  of  A.    This  system  of  equations  expresses  gross 
output  requirements  (x)  as  a  function  of  final  demand  (y)  and  the  tech- 
nological structure  of  the  economy  (A);  changes  in  the  elements  of  A  can 
come  about  for  a  variety  of  reasons.    In  this  paper  we  present  techniques 
for  solving  such  large  systems  of  equations,  and  for  updating  the  solution 
to  account  for  changes  in  A.   The  methods  presented  effect  substantial 
savings  in  computing  time  and  storage  requirements  over  those  convention- 
ally employed. 
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1.   INTRODUCTION 

In  an  input-output  econometric  model  we  are  often  concerned  with 
solving  the  system  of  n  equations  (I  -  A)  x  =  y,  repeatedly  for  various 
changes  in  the  elements  of  A.   This  system  of  equations  expresses  gross 
output  requirements  (x)  as  a  function  of  final  demand  (y)  and  the  tech- 
nological structure  of  the  economy,  and  changes  in  the  elements  of  A  can 
come  about  for  a  variety  of  reasons.    In  this  paper  we  present  techniques 
for  solving  such  large  systems  of  equations,  and  for  updating  the  solution 
to  account  for  changes  in  A.   The  methods  presented  effect  substantial 
savings  in  computing  time  and  storage  requirements  over  those  convention- 
ally employed. 

Several  methods  [l]  can  be  used  in  solving  the  system 

(I  -  A)  x  =  y  (1.1) 

1   3 
Gaussian  elimination  requires  —  n  multiplications,  Householder  triangular- 

2  3 
ization  requires  —  n  multiplications  and  n  square  roots,  and  Givens  tri- 
ll 3  12 
angular ization  requires  —  n  multiplications  and  —  n  square  roots.   In 

this  paper  we  use  a  modification  of  Givens  method  developed  by  Gentleman  [2] 

3 
that  requires  only  n  multiplications  and  no  square  roots.   There  are 

two  reasons  for  such  a  choice.   The  first  is  that  the  matrices  B  =  I  -  A 

are  usually  dense  and  of  large  size,  hence  on  some  computers  we  may  have 

to  resort  to  secondary  storage;  and  since  these  input-output  matrices 

are  often  stored  by  rows  Givens  transformations  are  more  suitable  than 

those  of  Householder.   Second,  some  of  the  procedures  in  updating  the 

solution  due  to  changes  in  B  such  as  removing  a  row  or  a  column  are 

inherently  unstable,  so  in  order  to  minimize  such  effects  it  is  desirable 

to  use  orthogonal  transformations. 


The  system  (l.l)  is  solved  by  Givens  transformations  as  follows. 
We  construct  an  orthogonal  matrix  Z  as  the  product  of  Givens  transforma- 
tions such  that  [2], 

ZB  =  D2R  (1.2) 

where  D  is  a  diagonal  matrix,  and  R  is  unit  upper  triangular.   From 
(l.l)  and  (1.2)  we  have 

D*Rx  =  Zy 
i.e.  , 

Rx  =  y  (1.3) 

which  is  solved  "by  backsubstitution. 

Note  that  we  do  not  store  Z,  the  Givens  transformations  are  multiplied 
by  the  final  demand  vector  y  as  they  are  produced.   When  the  matrix  B 
is  well-conditioned  instead  of  solving  (l.l)  we  may  solve  the  normal 
equation 

BtBx  =  Bty  (l.U) 

Of  course  the  condition  number  of  the  coefficient  matrix  is  the  square 
of  the  original  one;   however,  this  is  of  no  serious  consequence  since 
B  is  very  well-conditioned.   Using  (1.2)  we  may  write  (l.U)  in  the  form, 

RtDRx  =.  B*y  (1.5) 

Note  that  we  do  not  store  Z,  the  only  storage  required  is  that  for 

2 
B,  R,  D,  x,  and  y,  i.e.   2n  +  3n  words.   Finally  x  can  be  obtained  by 

solving  sequentially  the  systems, 

Rtx2  =  bV,   Dxx  =  x2,   Rx  =  x±  (1.6) 

When  a  new  right  hand  side  y  is  given  with  B  unchanged,  (1.6)  is 
again  used  where  we  replace  B  y  by  By. 

Various  methods  have  been  developed  for  updating  matrix  factorization 
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due  to  changes  in  the  elements  of  the  matrix.   We  refer  mainly  to  the 
papers  "by  Gill  and  Murray  [3]  and  Golub,  et.  al.  [h].      Throughout  this 
paper  we  use  Gentleman's  square  root  free  Givens  transformation  for  the 
updating  of  the  solution  x  in  (l.l)  due  to  row,  column,  and  a  single 
element  changes  in  B. 
2.   SQUARE  ROOT  FREE  GIVENS  REDUCTION 

We  describe  briefly  the  factorization  (1.2).   Let  us  consider  a 
Givens  transformation  [2]  that  rotates  two  row  vectors  such  that 


c   s 
-s  c 


va   u. 


va   uv 


/g  vx  /g"  v2  /a"  vn 

m  u 


y/oT       m  u 
0    /§"  v 


/s"  v 


v 


n_J 


(2.1) 


where 


a  =  a  +  3v 
3  =  a3/a 
c  =  a/a 


s  =  3v1/a 


u.  =  cu.  +  sv. 
Ill 


V.  =  V.  -  V.     u. 

1    1    1  1 


i  =  2,.  .  . ,  n 


(2.2) 


Note  that  no  square  root,  evaluations  are  involved  in  these  formulae  and 
(2.1)  may  he  written  as 
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We  define  an  n-dimensional  Givens  transformation  Z   by, 

K. 


•i 


-s 


*1 


i-th  row 


k-th  row 


J 


(2.1t) 


in  which  c  =  cos  a  . ,  s  =  sin  a  . ,  and  the  angle  a  .  is  chosen  such  that 
the  element  in  the  position  (k,  i)  of  the  matrix  Z,   B  is  eliminated.   Let 


2  =  *li 


z1  ,  z1 

n-1  n 


and. 


thus , 


Z  =  Z 


n-1 


k 


Z2Z1 


(2.5) 


ZB  =  D"R 

where  D  is  diagonal,  and  R  is  a  unit  upper  triangular  matrix. 

2  3 

As  is  mentioned  in  Section  1,  2n  +  3n  storage  locations  and  n 

multiplications  without  square  root  operations  are  necessary  for  Givens 

reduction.    In  this  paper  and  the  attached  programs  we  are  assuming  that 

the  matrices  B  and  R,  and  the  diagonal  vector  D  are  stored  in  the  high  speed 

memory  of  the  computer  and  are  available  for  updating  the  solutions  due  to 

changes  in  B,  which  will  be  discussed  in  the  following  sections. 

3.   COLUMN  MODIFICATION 

Given  the  system  (l.l)  and  the  factorization  (1.2)  we  wish  to  determine 

the  new  factorization 

B^B  =  R^DR  (3.1) 
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in  less  than  n  operations,  where  B  is  different  from  B  in  only  one  column 

([31,  [W. 

We  first  consider  deleting  a  column  b.  from  B.   Let 

\°   [V-"'bi-l'  bi+l«  "•'  V 

be  an  n  x  (n-l)  matrix  obtained  from  B  by  deleting  the  i-th  column,  thus 

t       t 
B,   B.  =  R  D  R 
1   1    o    o 

where  R  is  R  with  the  i-th  column  deleted.    If  we  partition  R  as 
o  *  o 


R  = 
o 


H, 


H, 


then  the  (i-l)  x  (i-l)  submatrix  H  is  unit  upper  triangular,  and  H  is 
an  upper  Hessenberg  matrix  of  order  (n-i).  We  wish  to  find  an  orthogonal 


transformation  Z  such  that 


B*  B,  =  B}   dVzD^R 
1   1    o        o 

"  EiDiRi 


(3.3) 


where 


k 


k 


D,  %  =  ZITR 
1  1       c 


in  which  R  is  an  (n-l)  x  (n-l)  unit  upper  triangular  matrix.   This  c 


an 


be  done  by  defining  Z  as  a  product  of  the  following  Givens  transformations: 

z  =  z11"1.  .  .  .  zni+i  z1  n 

n  1+2   i+1 

That  is,  Z  reduces  the  upper  Hessenberg  matrix  H  to  a  unit  upper  triangu- 
lar matrix.    Thus,  deleting  a  column  from  B  and  obtaining  a  modified 
factorization  involve  only  the  temporary  storage  of  the  (n-i)  x  (n-i)  lower 
principal  submatrix  of  R  . 
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We  now  consider  adding  a  column  b  to  B  . 


Let 


B  = 

[V 

*] 

R   = 

A 

r 

0 

_ 

(3.U) 


and 


D  = 


in  which  an  n-vector  r  with  its  n-th  component  unity  and  a  scalar  d  are 


to  be  determined  such  that 


B^B  =  RtDR 


From   (3.U), 


B*B- 


A 


\\ 


b*fc 


B      .'      b 

i  : 


V 


b*b 


(3.5) 


and 


eH  R  = 


~D1 

—1 

:    o 

~Ri 

r 

0 

:    d  _, 

_   0 

w 

[h\ 


0    ]   r 


i^Dr 


(3.6) 


Comparing  right  hand  sides  of  (3-5)  and  (3.6),  we  obtain 
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r  =  B  b 


(3.7) 


r  Dr   =  b  b 
from  which  r  and  d  are  computed. 

Thus,  the  modified  factorization  (3.1)  is  completely  determined  and 


the  modified  solution  x  is  obtained  from 

RtD  R  x  =  bV 
by  forward  and  backward  substitution  as  before 


(3.8) 


The  maximum  number  of  operations  necessary  for  a  column  modification 

2 
is  only  5n  multiplications.   Note  that  only  3n  extra  storages  for 

x,  y,  and  a  new  column  in  addition  to  those  for  B,  D,  and  R  are  required 

for  deleting  and  adding  a  column,  and  updating  D  and  R. 

k.      ROW  MODIFICATION 

Let  us  first  consider  deleting  a  row  b.   from  B, 

l        ' 


Bi  =  B- 


then 

B,\    =  B*B  -  b.b* 
11         11 

=  R*DR  -  b.b* 
l  i 

where  R  and  D  are  the  original  factors  of  B  and  we  assume  that  they  are 

available.   We  wish  to  determine  a  unit  upper  triangular  matrix  R.  and  a 


diagonal  matrix  D  such  that 


R^D-R-  =  R^R-b.b? 
Ill        li 


To  do  this  ([3],  [U]),  solve 


R  Dv  =  b, 


(U.l) 


(U.2) 


for  v  and  let 
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2        i  II  2 

s     =  1   -  v     2 

Nov  we  construct  an   (n+l)   x   (n+l)  matrix 


R     = 
o 


D2v 


D%R 


-  1 

0~] 

s 

0™ 

0 

_ 

V 

R 

and  reduce  it  to  a  unit  upper  triangular  matrix, 


"- 

1 
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:    ° 
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0 

_  0 

■      % 
D 

V 
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0 

1 

o    \ 
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1  -J 


(fc.3) 


where 


Z  =  z: 


Z1   Z1.. 
n   n+l 


and,  R  and  D  are  the  desired  factorization  of  B  .    In  order  to  show 


that,  we  have  from  (k.3. 


2  u  tn 
s   +  v  Dv 


vtDR 


RtDv 


r'Sr 


1 

t 
r 

r 

rr*  +  R1tD1R1  ^ 

comparing  both  sides,  we  obtain 


vtDv  +  s2  =  1 


RtDv 


=  r 


Therefore 


R^R      =  rr*  +  R-^DjR-, 


r  =  b. 

l 


and,  the  relation  (U.l)  is  satisfied. 


For  adding  a  row  g.  to  B  ,  let 


B  =  B  + 
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Then 


B^B  =  B1tB1  +  g.g^ 


(k.k) 


Now  construct  an  (n  +  l)  x  n  matrix 


\\ 

t 

.  gi  J 

pper  triangular  matri 

DiRi 

D'2R 

t 
gi 

0 

(U.5) 


where 


From   {k.5)  we  obtain 


Z  =  Z\i       '    '    '      z2+i      Z\n 
n+1  n+1       n+1 


Rl\Rl  +  gigit   =  *H* 


(k.6) 


Therefore,  from  (k.k)   and  (k.6) 


£*£  =   ^DR 
Thus  the  modified  factorization  is  completely  determined  and  the  modified 
solution  is  obtained  from 

R^DRx  =  bV 
as  before. 

The  number  of  operations  required  for  the  row  modification  is 

2 
5%n  multiplications.   Again,  only  Un  extra  storage  locations  are 

necessary  in  addition  to  those  for  B,  D,  and  R. 

For  an  alternative  method  for  row  modification  see  Sameh  and 

Bezdek  [5]. 
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5.   ELEMENT  MODIFICATION 

When  only  one  element  of  the  matrix  B  is  changed  and  we  wish  to 

update  the  solution,  the  Sherman-Morrison  formula  [6]  may  he  utilized 

2 
since  only  n  multiplications  are  required  in  the  process. 

Let  us  consider  the  prohlem  in  which  we  wish  to  solve 

(B  +  a  e  e  )  x  =  y  (5.1) 

where  a  is  a  scalar,  and  e.  and  e.  are  the  i-th  and  the  j-th  column 

1  -J 

vectors  of  an  identity  matrix,  respectively. 

From  (5.1)  and  the  Sherman-Morrison  formula,  we  obtain 

x  =  (B  +  ae.  e .  )"  y 

1   J 

=  (B_1  -  SB-1  e.  e.1  B_1 )  y 

=  B_1y  -  -3  B"1  e  e  *  B_1y  (5-2) 

where  3  =  (a_1  +  e  }   B_1  e.  )  _1  (5-3) 

j       1 

Assuming  the  original  solution  is  known,  x  =  B  y,  (5.2)  becomes 

x=x-3B"1e.e.tx  (5.U) 

1   j 


If  we  solve  for  w  in 


B  w  =  e±  (5.5) 


using  the  factorization  (l.5)»  then 

3  =  (a-1  +  w. )-1 
J 


and 

X=x-3x.w  (5-6) 

J 

2 
Although  n  multiplications  are  necessary  for  the  element  modifications, 

if  the  i-th  column  of  B   is  available,  the  number  of  multiplications  is 

reduced  to  n. 

Only  2n  extra  storage  locations  are  necessary  in  addition  to  those 

for  B,  D,  and  R. 
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Until  now,  we  have  as  sinned  that  the  computer  core  memory  is  large 
enough  so  that  the  matrices  B,  R,  and  D  can  simultaneously  be  contained 
in  it.   However,  since  we  are  dealing  with  matrices  of  large  size,  it 
may  happen  that  the  core  memory  cannot  accommodate  all  or  part  of  them 
simultaneously.   In  this  case,  the  matrices  can  be  stored  by  rows  in 
the  secondary  storage,  and  several  rows  of  B  or  R  are  brought  into  the 
core  whenever  necessary. 

For  applying  Givens  transformation  for  triangularization  of  B  or 
updating  R,  only  as  many  rows  of  B  or  R  that  the  core  memory  can 
accommodate  are  brought  in  since  only  two  rows  of  B  or  R  are  necessary 
for  a  single  Givens  transformation. 

In  order  to  compute  B  y,  which  appears  in  (3.7),  (3.8),  and  Section  k9 
each  row  vector  of  the  matrix,  B.  in  the  core  is  multiplied  by  y. ,  i-th 


component  of  y,  i.e. 


t    n 
B  y  =  Z  B.y 

i=l 
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SUBRGUTINt  QIvQEN<mNj  N,  A,  B,  C>,  R..  Y,.<  EPS> 
C 

C     THIS  SUBROUTINE  SOLVES  fi  SVSTEM  OF  LINEAR  EQUATIONS  A*X=B  BY 
C     GIVENS  TRANSFORMATIONS  WITHOUT  SQUARE  ROOT  EVALUATIONS. 
C      IT  DETERMINES  THE  FACTORIZATION  Z*B=D.**<1/2>*R,  WHERE  Z    IS  AN 
C     ORTHOGONAL  MATRIX,  D  IS  A  DIAGONAL  MATRIX..  AND  R  IS  A  UNIT  UPPER 
C     TRIANGULAR  MATRIX.  Z  IS  NOT  STORED   SOLUTION  IS  THEN  OBTAINED 
C     BV  A  BACKWARD  SUBSTITUTION.  IF  ANY  ELEMENT  OF  A  IS  LESS  THAN  EPS> 
C     IT  IS  REPLACED  BV  ZERO. 
-C 

C     INPUT 

C        NN   DECLARED  DIMENSION  OF  MATRIX  A 
C        N    ORDER  OF  MATRIX  A 
C        A    MATRIX  OF  ORDER  N 
C        B     RIGHT  HAND  SIDE  VECTOR 
C 

C     OUTPUT 

C        X    SOLUTION  VECTOR 
C        D     DIAGONAL  VECTOR 

C        R    UNIT  UPPER  TRIANGULAR  MATRIX  OF  ORDER  N 
C        EPS   NACEPS+CNORM  OF  A> 
C 

P.EAL+3  A<:NN,  NX.  B<N>,  D<N>,  RCNN,  N>,  X<N> 

REAL+3  ALPHA,  BETA,  UI,  VI,  CBAR,  SEAR,  TEMP,  ANORM,  MACHEP,  DMAX1, 
X       DABS, EPS 

DATA  MACHEP/Z3410000000000000/' 

DO  5  J=l,  N 

DCJ>=1.  0D9 

X<J)-BCJ5 

DO  5  1=1,  N 
5  R(LJ)=fl(I,J) 
C 
C     COMPUTE  INFINITE  NORM  OF  A  AND  EPS=MACHEPS+ ANORM 

RNORM=0.  OD0 

DO  15  1=1,  N 

TEMF-O.  0D0 

DO  10  J=l, N 
10  TEMF-TEMP+DABSCRC 1, J) > 
15  ANORN=DMAXl<: ANORM,  TEMP) 

EPS=MACHEF+ANORM 
C 

NM1=N-1 

DO  180  J=l,  NM1 

JP1=J+1 

NJ=N+JP1 

I=J 

UI=R<I,J) 
C 
C     SPECIAL  CASES  FpR  GIVENS  TRANSFORMATIONS 

IF  <DABS<IJI>.  LT.  EPS)  GO  TO  25 

DO  20  H=.JP1,  N 
28  R(I,  M>=R'r.I,M)/UI 

X<I)=X<IVUI 

R(I,  J)=l.  OD0 

DCI)=UI*UI*D<I> 

GO  TO  38 
25  R<I,  J)=0.  0D0 
38  CONTINUE 

DO  100  L=JP1, N 

K=NJ-L 
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V1=K',K.-  J.> 

IF    (DABS^'VI).  GE.  EPS)    GO   TO   35 

ROG  J)  =8.  6D0 

GO   TO   39 
25    CONTINUE 

IF    CDRBSail).  GE.  EPS)    GO   TO   58 

DO   40   M--J,  N 

TENP=-F:a,M> 

Re:  i,  m)=r<k,  M) 
48   R<K,  M)=TEHP 

TEMP=-X(I) 

Xa>=X<K) 

X<K)=TEMP 

TEMP=D<I) 

D<I)=D<K) 

DOO=TEMP 

DO   45   M=JP1,  N 
45   R<CI,M)=Ra,M)/VI 

X(I)=X(I 5/VI 

DCI)=VI+VI+D<I) 

R(L  J)=l.  8D8 

GO  TO  58 
58  CONTINUE 
C 
C     GIVENS  TRANSFORMATIONS 

ALPHA=D<I) 

beta=dc:k> 

d< i >=alpha+beta+vi+vi 

DaO=RLPHA+BETfl/D<:  I ) 

R(L,T)=1  0D8 

R(K,  J)  =8.  8D3 

CBAR=RLPHA/D<I) 

SBRR=BETfl*VI/D<I) 

DO   55   M=JP1,  N 

TEMF-RaC  M)-VI*R<L  M) 

R<  I,  N::'=CBAR*R<:i,  M>+SBflR+R<K,  M) 
55  RCK,  r-D=TEMP 

TENP=HOO-VI+Xa) 

Y,(.  I  >=CBRR+Xa  >+SBAR+X<K) 

X<K)=TEMP 
98   CONTINUE 
188   CONTINUE 

TEMP=R<N, N) 

X<N>=X'::N>/TEMP 

D<N)=TEMP*TENP*D<N) 

RCN,  H)=l.  8D8 
C 
C     BACKWARD  SUBSTITUTION 

DO  128  1=1,  N 

J=N-I+1 

JP1=J+1 

TEriP=8.  8D8 

IF  CJP1.  GT.  H)  GO  TO  128 

DO  118  n=JPl, N 

118  tenf-temp+rcj,  m)*x<m) 
128  k<j>=x<j:>-tenp 
138  continue 

RETURN 

END 

SUBROUTINE  COLNOD<NN,  N,  A,  B,  COL,  ICOL,  D,  R,  X,  EPS) 
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THIS  SUBROUTINE  UPDATES  THE  SOLUTION  X,  DIAGONAL  MATRIX  D,  AND 
UNIT  UPPER  TRIANGULAR  MATRIX  R  NHEN  A  COLUMN  OF  A  IS  CHANGED, 
PROVIDED  THAT  THE  FACTORS  P.  AND  D,  OBTAINED  FROM  'GIVGEN', 
ARE  KNOWN. 


INPUT 
NN 
N 
A 
B 

COL 
I  COL 
D 
R 
EPS 

OUTPUT 
v 

D 
R 


DECLARED  DIMENSION  OF  MATRIX  A 

ORDER  OF  MATRIX  A 

MATRIX  OF  ORDER  N 

RIGHT  HAND  SIDE  VECTOR 

NEW  COLUMN  VECTOR 

THE  ICOL-TH  COLUMN  OF  A  IS  CHANGED 

DIAGONAL  VECTOR  OBTAINED  FROM  ■•'GIVGEN' 

UNIT  UPPER  TRIANGULAR  MATRIX  OBTAINED  FROM  'GIVGEN' 

MACHEPS+CNORH  OF  A';-  OBTAINED  FROM  ''GIVGEN' 


UPC 
UPC 
UPC 


ATED  SOLUTION  VECTOR 
ATED  DIAGONAL  VECTOR 
ATED  UNIT  UPPER  TRIANGULAR  MATRIX 


RE AL+8  A <  NN,  N > ,  B < N > ,  COL < N ) ,  D < N > ,  R < NN,  N > ..  X < N > 
REAL+S  ALPHA,  BETA,  UI,  VI,  CBAR,  SBAR,  TEMP,  EPS,  DABS 

REMOVE  ICOL-TH  COLUMN  OF  R  TO  OBTAIN  R9 

NM1=N-1 

IF  a  COL.  EQ.  W    GO  TO  218 

DO  20  J=ICOL, NN1 

JP1=J+1 

DO  20  1=1,  JP1 

RU,  J)=R<L  JP1) 


TRIANGULAR  FORM  Rl 


REDUCTION   OF   RO   TO   UPPE 

DO   208    I=ICOL, NM1 

K=I+1 

UI=R<I, I) 

VI=R<K, I> 

IF    CDABSailX  LT.  EPS)    GO 

IF    (I.  EG.  NMD    GO   TO   115 

DO   110   J=K, NM1 

lis 

R<I,  j>=r<i..  jvui 

115 

RCI,  I)=l.  OD0 

Da)=UI*UI*DCI) 

GO   TO   125 

120 

R(I,  I>=8.  8D0 

125 

CONTINUE 

IF    ':DABS<VI>.  GE.  EPS)    GO 

R<K,  I>=0.  0D0- 

GO   TO   190 

130 

CONTINUE 

IF    CDABSaJIX  GE.  EPS>    GO 

DO   135   J=I,NM1 

TEMP--R<I, J) 

R<I,  J>=R<K,  J) 

135 

ROC  J>=TEMP 

TEMF-DCI) 

D(I>=D<K) 

DOO=TEMP 

IF    (I.  EQ.  NM1>    GO   TO   145 

TO  128 


TO  139 


TO  150 
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DO  i4G  J=K..  Nr-ii 

140  R  <.  I ..  J 1  =R  <  I ,  ,T  )  / V I 
145  CONTINUE 

DCI)=VI*VI*D<I) 

R(L  I>=1.  ODO 

GO  TO  ISO 
150  CONTINUE 

RLPHR=D<I> 

BETfl=D<K> 

D  C I  :<  =RLPHR+BETR+V  I +V I 

D  i.  K  )  =RLFHR  +  EETfl/D  <  I  > 

R(I..  I)=l.  ODO 

RCK,  I)=0.  ODO 

CBRR=flLPHfl/DCO 

bBRR=BETH+VI/Da> 

IF  a.  EQ.  NM1>  GO  TO  190 

DO  155  .J=K,  NM1 

TEMP=R < K,  J  > - V I +R  C I ,  J ) 

R<  I,  J>=CBRR*R<  I,  J>+SBRR*P<K.  J) 
155  RCK,  J."J=TEMP 
190  CONTINUE 
200  CONTINUE 
C 

C      CTRRNSPOSE  OF  fll)*COL=X 
210  CONTINUE 

IC0LM1=IC0L-1 

IF  (I COL.  EQ.  1>  GO  TO  235 

DO  230  J=l,  IC0LM1 

TEMP=0.  ODO 

DO  220  I  =1,  N 
220  TEMP=TEMP+Ra,  J>+COL<  I  ) 
230  K<J)=TEMP 

IF  a  COL.  EQ.  N>  GO  TO  260 
235  CONTINUE 

DO  250  ,T=ICOL,  NM1 

•JP1=J+1 

TEMP=0.  ODO 

DO  240  1=1,  N 
240  TEMP=TEMP+R<I,  JP1 > +COL < I> 
250  K<:J>=TEMP 
260  CONTINUE 
C 

COMPUTE  THE  LRST  COLUMN  OF  RTILDE  AND  THE  LRST  COMPONENT  OF  DTILDE 

DO  320  J=l,  NM1 

JM1=J-1 

TEMF-O.  ODO 

IF  <J.  EQ.  1)  GO  TO  320 

DO  310  1=1, JM1 
310  TEMP=TEMP+Ra,  J)+R<  I,  N) 
320  PC  J, N)=X<J)-TEMP 

DO  330  1=1,  NM1 
330  R<I,N)=R<I,N)/D<I) 

TEMP=0.  ODO 

DO  340  1=1,  N 
340  TEMP=TEMP+COL< I >*COL< I ) 

BETfl=0.  ODO 

DO  350  1=1,  NM1 
350  BETR=BETR+R< I,  N  VtRC I, N)#D< I > 

D<N)=TEMP-BETfl 

R<N,  N>=1.  ODO 
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u 

C     COMPUTE  RHS=<RTILDE)T+B 

IF  <I COL.  EQ.  1>  GO  TO  235 

DO  380  J=l,  IC0LM1 

TEMP=9.  ODO 

DO   376    1=1..  N 
370   TEMF-TEMP+FKI,  J>+B<I) 
380  :>;•::  J  >=TEMP 

IF  (I COL.  EQ.  N)  GO  TO  405 
385  CONTINUE 

DO  400  J=ICOL,  NM1 

JP1=J+1 

TEMP=0.  ODO 

DO  330  1=1,  N 
390  TEMP=TEMP+fi< I,  JF1>*B< I > 
40O  X<J)=TEMP 
405  TEMF-O.  ODO 

DO  410  1=1,  N 
410  TENF -TEMP+COL C I ) +B  < I ) 

X<M>=TEMP 
C 
C     FORWARD  RND  BACKWARD  SUBSTITUTIONS 

DO  430  J=l,  N 

JM1=J-1 

TEMF-O.  ODO 

IF  CJ.  EQ.  1>  GO  TO  430 

DO  420  1=1,  JM1 
420  TEMP=TEMP+R< I,  J>*X< I > 
430  X<J>=X<J)-TEMP 

DO  440  1=1,  N 
440  X(I)=X(IVD(I) 

DO  460  1=1, N 

J=N-I+1 

JP1=J+1 

TEMF=0.  ODO 

IF  CJP1.  GT.  N)  GO  TO  460 

DO  450  K=.JP1,  N 

450  temf-temp+rcj,  k>+x<k> 
460  xcj:>=xc:j;>-temp 

C 

C     REORDERING  OF  Xtl), .  .  .  ,l«H> 

IF    U COL.  EQ.  N>    GO   TO   480 

TENF-XaD 

DO   470    1  =  1 COL,  NM1 

J=N-I+ICOL 
470  X<J)=X<J-1> 

XaCOL>=TEMP 
480   CONTINUE 

RETURN 

END 

SUBROUTINE  ROWMOD<NN,  N,  R,  B,  ROW,  IRON,  D,  R,  X,  EPS,  V) 
C 

C     THIS  SUBROUTINE  UPDATES  THE  SOLUTION  X,  DIFlGONflL  MATRIX  D,  RND 
C     UNIT  UPPER  TRIANGULAR  MATRIX  R  WHEN  R  ROW  OF  MATRIX  R  IS  CHRNGED, 
C     PROVIDED  THRT  THE  FACTORS  R  RND  D,  OBTAINED  FROM  'GIVGEN', 
C     RRE  KNOWN.  V  IS  A  TEMPORRRV  STORAGE  VECTOR. 
C 

C     INPUT 

C        NN   DECLARED  DIMENSION  OF  MATRIX  fl 
C        N    ORDER  OF  MATRIX  fl 


-IT- 


'-:  h  fiHiRIH  OF  ORDER  N 

C        B     RIGHT  HAND  SIDE  VECTOR 

C         RON   MEN  ROW  VECTOR 

C         I ROW  THE  IROW-TH  ROW  OF  fl  IS  CHANGED 

D     DIAGONAL  VECTOR  OBTAINED  FROM  'GIVGEN' 

UNIT  UPPER  TRIANGULAR  MATRIX  OBTAINED  FROM  'GIVGEN' 

EPS   MACHEPS+<NORM  OF  A>  OBTAINED  FROM  'GIVGEN' 

C     OUTPUT 

C        X     UPDATED  SOLUTION  VECTOR 

C        D     UPDATED  DIAGONAL  VECTOR 

R     UPDATED  UNIT  UPPER  TRIANGULAR  MATRIX 
C- 

REAL+3  A  <:  NN,  N  >  >  B  <  N  ) ,  ROW  <  N  ) ,  D  <  N  ) ,  P  '■  NN  •  N > .  X < N ) ,  V ' N> 

REAL+8  ALPHA,  BETA,  UI,  VI,  CBAR,  SBAR,  TEMP,  DI,  DK,  S,  EPS,  DABS 

C     SOLVE  <RT>*<CD>+V=BI  FOP  V 
C 

DO  10  J=l,  N 
10  XOJ)=A<IR0W,  J) 
DO  30  J=l,  N 
JM1=J-1 
TEMF-O.  OD0 

IF    <..].  EQ.  1)    GO    TO   30 
DO   20    1=1,  .TNI 
20   TEMP=TEMP+Ra,  J)#V<  I  ) 
38    V<J>=XCJ>-TEMP 

DO   40    1=1,  N 
40   V(I)=V(DVD(I) 
C 

C     REDUCTION  OF  R0  TO  UNIT  UPPER  TRANGULAP  MATRIX  PI 
DO  45  1=1,  N 
45  X(  I  )=8.  0D9 
DI=1.  OD0 

S  =  0.    0D0 

DO   20O    1=1,  N 

K=N-I+1 

UI=S 

VI=VOO 

IF  <DflBS<VIJ.  GE.  EPS)  GO  TO  135 

VaO=0.  0DO 

GO  TO  130 
135  CONTINUE 

IF    (DABSCUi:.'.  GE.  EPS)    GO   TO   150 

DO   140   J=K, N 

TEMP=-X<J> 

X<J>=RCK,  J) 
140   ROG  J)=TEMP 

TEriP=DI 

DI=D<K> 

DOO=TEMP 

DO   145   J=K,  N 
145   X(.J)=XCJ)/VI 

DI=VI*VI*DI 

5=1.  OD0 

V<K)=9.  OD0 

GO  TO  190 
150  CONTINUE 

ALPHA-DI 

BETA=DOO 
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L>  1  ""HLPHN+bh  !  H+V1+V1 

D  <  K >=flLPHfl+BETfl/D I 

5=1.  0D0 

V(K>=8.  QD0 

CBAR=ALPHA/DI 

SBAR=BETA*VI/DI 

DO  160  J=K..  N 

TEHF-ROG  J>-VI*X<J> 

X  '■  -J  )  =CBAR#X  <.  J  >  +SBAR*R<K,  J  ) 
160   RQG  ,J>=TEMP 
190   CONTINUE 
2O0   CONTINUE 
C" 
C-  COTiF-UTE  RH5=  CAT  ILDEIm^B 

DO   220   J=l, N 

TEMF-0.  OD0 

DO  210    1=1,  N 

S=A<I,  J) 

IF  <I.  EQ.  IROW>  S=ROW<J> 
210  TEMP=TEMP+S+B< I > 
220  X<J>=TEMP 


C 


"¥2 


REDUCTION  OF  R2  TO  UNIT  UPPER  TRIANGULAR  MATRIX  RTILDE 

DK=1.  OD0 

DO  400  1=1,  N 

IP1=I+1 

UI=R<I, I> 

IF  CDABSOJIXLT.EPS)  GO  TO  325 

IF  (I.  EQ.  N)  GO  TO  322 

DO  320  J=IP1, N 

R(LJ)=R(I,  J)  All 

CONTINUE 

R<I,  I>=1.  OD0 

D<i>=UI*UI*DCI) 

GO  TO  330 
325  R<I,  I>=9.  0DO 
330  CONTINUE 

VI=ROH<I) 

IF  C  DABS  <  VI).  GE.  EPS)  GO  TO  335 

ROW  CI  5=0.  OD0 

GO  TO  390 
335  CONTINUE 

IF  <DABSaU>.  GE.  EPS)  GO  TO  350 

IF  <I.  EQ.  N)  GO  TO  342 

DO  340  J=IP1,  N 

TEMP=-R<I,  J) 

R<I,  J)=ROl-KJ) 
340  ROW<J)=TEMP 
342  CONTINUE 

TEf1P=D<I) 

DU>=DK 

DK=TEMP 

if  a.  eq.  n:>  go  to  347 

DO  345  J=IP1, N 
345  RCI,  J)=R(LJ)/VI 
347  CONTINUE 

D<I)=VI*VI*D<I> 

R<I,  I)=l.  0D0 

GO  TO  390 
350  CONTINUE 
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c 


RLPHR=D<I> 
BETfl=DK 

DC  I >=RLPHR+BETR+VI+VI 

DI^ALFHA+BETA/Dd) 

CBRR=RLFHR/D<i:> 

SBRR=BETR+VI/DCI) 

R(L  I  ;<=!.  0D8 

ROW<n=0.  0D0 

IF  H.  EQ.  N)  GO  TO  3'30 

DO  360  J=IP1,  N 

TENF -ROW C J  > - V I *R ( I ,  J ) 

R  C I ,  J  )  =CBP,R+R  (LJ  :>  +SBRR+ROW  <  J  ) 
360  RONCJ)=TENP 
330  CONTINUE 
400  CONTINUE 

FORWARD  FIND  BRCKWRRD  SUBSTITUTION 
DO  430  J=l,  N 
JM1=J-1 
TENF-O.  OD0 

IF  CJM1.  LT.  1)  GO  TO  430 
DO  420  1=1,  ..TNI 
420  TEMP=TEMP+R  <  I ,  J  ;< *V  <  I  > 

430  v<:.n=v  0:0 -temp 

DO  440  1=1,  N 
440  VCI)3V(I)/D<I) 

DO  460  1=1,  N 

J=N-I+1 

JP1=.J+1 

TENP-0.  0D0 

IF  <JP1.  GT.  ID    GO  TO  460 

DO  450  K=.JP1,  N 
450  TEMF -TEHP+R < J,  K > +X < K ) 
460  X<J>=VCJ)-TEMP 

RETURN 

END 

SUBROUTINE  EMTMODCNN,  N,  fl,  X,  D,  R,  EMT,  ITH.  JTH,  V> 


C  THIS  SUBROUTINE  UPDATES  THE  SOLUTION  X  OF  fi+X=B  WHEN  RN  ELEMENT 

C  OF  R  lb  CHANGED,  PROVIDED  THAT  THE  FACTORS  R,  D  AND  THE  SOLUTION 

X,  OBTAINED  FROM  "  G IVGEN' ,  ARE  KNOWN.  V  IS  A  TEMPORRPV  STORAGE 

C  VECTOR. 
C 

C  INPUT 

C  NN    DECLARED  DIMENSION  OF  MATRIX  R 

C  N     ORDER  OF  MATRIX  fi 

C  A     MATRIX  OF  ORDER  N 

C  X    SOLUTION  VECTOR  OBTAINED  FROM  'GIVGEN' 

C  D     DIAGONAL  VECTOR  OBTAINED  FROM  'GIVGEN'' 

C  R     UNIT  UPPER  TRIANGULAR  MATRIX  OBTAINED  FROM  'GIVGEN' 

C  EMT   NEW  ELEMENT 

c  ITH, JTH   THE  < ITH,  JTH > -ELEMENT  OF  R  IS  CHANGED 
C 

C  OUTPUT 

C  X    UPDATED  SOLUTION  VEPTAR 

C 

REAL+3  A  CNN,  N),  X(N>,  D<N),  RCNN,  NX.  V(N),  SIGMA,  TAU,  EMT,  TEMP 

C  SOLVE  RT+D+R+V=RT*E< I ) 

DO  3G  J=l,  N 
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TEMP=0.  8D0 

IF  <J.  EQ.  1)  GO  TO  30 

DO  20  1=1,  JM1 
20  TENP-TEMP+RCI,  JD+VC  I  ) 
30  V<J)=fiCITH,  J) -TEMP 

DO  40  1=1,  N 
40  V(I)=V(I)/D(I) 

DO  60  1=1, N 

J=N-I+1 

JP1=J+1 

TEMP=0.  0DO 

IF  <J.  EQ.  N)  GO  TO  60 

DO  50  K=JP1, N 
50  TEMF-TEMP+ROJ,  K>*V<K> 
60  V<J)=V<J>-TEMP 
C 
C     CONFUTE  SIGMA  RND  TRU 

SIGMfl=EMT-naTH,  JThD 

TflU=l.  0DO+S I GNR+V  < JTH  > 

TRU=SIGMR/TRU 
C 
C     COMPUTE  THE  UPDATED  SOLUTION 

TEMF-TRU+X'r...TTH> 

DO  70  1=1,  N 
70  X  <- 1  >  =X  <  I  >  -TEMP*V  <  I ) 

RETURN 

END 
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